# Stratified, multilevel sampling
# J. Chen & M. Walsh (Jul. 2012)

# Required packages, to install use:
# install.packages(c("sampling","rgdal","proj4","lattice"), dependencies=TRUE)
require(sampling)
require(rgdal)
require(proj4)
require(lattice)

# Load data
# setwd("~/Dropbox/ethiopia/Priority Woredas/sampling_code")
sampling_centers <- read.csv("Priority_Ag.csv")
sampling_centers_CNS <- sampling_centers #[sampling_centers$Survey.by=="CNS", ]

#####################################
# inputs needed
# 1. total number of sampling points (N)
# 2. number of 100 m grid cells within each 500 m grid cell (N.with500m)
# 3. stratification variable (strata)
#####################################

N <- 4000 
N.with500m <-  2 
N.each1k <- N.with500m*4
strata <- "Woreda"

stratasize <- table(sampling_centers_CNS[, strata])
stratasize <- stratasize[stratasize>0]

#number of samples for each woreda
N.strata.1k <- floor((N/(N.each1k *sum(stratasize))*stratasize+0.5))

#set random seeds
seed.selected <- 20120710
set.seed(seed.selected)
sampled_locations <- rep(NA, 9)
strata.names <- names(stratasize)

# drawing a hierachical sample
for(i in 1:length(strata.names)){
	woreda.sample <- sampling_centers_CNS[sampling_centers_CNS[, strata]==strata.names[i], ]
	
	# sample the locations of centers of site center points
	if(N.strata.1k[i]>0){
		sentinel.sample <- sample(1:dim(woreda.sample)[1], N.strata.1k[i])
	
		# hierachical sampling for each site
		for(j in 1:length(sentinel.sample)){
			xcenter <- woreda.sample[sentinel.sample[j], "Easting"]
			ycenter <- woreda.sample[sentinel.sample[j], "Northing"]
			xoff <- xcenter -250
			yoff <- ycenter -250
			
		# xdim & ydim specifiy the number of cells to be sampled in x & y directions
			xdim <- 10 
			ydim <- 10

		# specify grid resolution (grain, in m)
			grain <- 100

	    # generate the grid
			grid <- as.data.frame(coordinates(GridTopology(c(xoff,yoff), c(grain,grain), c(xdim,ydim))))
			colnames(grid) <- c("x", "y")

		# set up level ID's at the desired scales (res.pixel, in m)
	  		res.pixel <- c(500, 100, grain) 
	  		X2 <- ceiling((grid$x-xoff+1/2*xdim)/res.pixel[1])
	  		Y2 <- ceiling((grid$y-yoff+1/2*ydim)/res.pixel[1])
			L2 <- cleanstrata(paste(X2, Y2, sep=""))

			X1 <- ceiling((X2*res.pixel[1]-(grid$x-xoff))/res.pixel[2])
			Y1 <- ceiling((Y2*res.pixel[1]-(grid$y-yoff))/res.pixel[2])
			L1 <- cleanstrata(paste(X1, Y1, sep=""))

			X0 <- ceiling((X1*res.pixel[2]+(grid$x-xoff-X2*res.pixel[1]))/res.pixel[3])
			Y0 <- ceiling((Y1*res.pixel[2]+(grid$y-yoff-Y2*res.pixel[1]))/res.pixel[3])
			L0 <- cleanstrata(paste(X0, Y0, sep=""))

		# update grid with level ID's
			grid <- cbind(grid, L2, L1, L0)

		# draw a sample
		# s.size specifies the sample size at each level (e.g. L2=4, L1=2, L0=1)
			s.size <- list(4, rep(N.with500m, 4), rep(1, N.with500m*4))
			s <- mstage(grid, stage=list("cluster", "cluster", ""), varnames=list("L2", "L1", "L0"), size=s.size, method="srswor")
			sample.DGG <- getdata(grid, s[[3]])
			plot1 <- xyplot(y~x, data=sample.DGG, xlab="Easting (m)", ylab="Northing (m)", pch=3, cex=0.5, asp=1)
			sampled_locations <- rbind(sampled_locations, cbind(sample.DGG, Woreda=strata.names[i]))
		}
	}
	print(i)	
}
sampled_locations <- sampled_locations[-1, ]

# reproject to Lat/Lon, WGS84 coordinate reference system
latlon <- project(cbind(sampled_locations$x, sampled_locations$y),  "+proj=utm +ellps=WGS84 +zone=37 +north=T +units=m +no_defs", inv=TRUE)
sampled_locations$x <- latlon[,1]
sampled_locations$y <- latlon[,2]

setwd("sampled_locations")
write.csv(sampled_locations, "samplelocations.csv")

#generate kml file
system("python spatial_csv_to_kml.py")
